function [W]= melfilter(fs, N, nofChannels)
%function [W,Ws]= melfilter(fs, N, nofChannels)
%calculo de banco de filtros Mel
%Paula Sanz Octubre de 2008


%parametros
%-----------------------------------------
%fs          = 16000;    %sampling frequency
%N           = 1024;     %FFT Window length
%nofChannels = 50;       %number of mel channels
%-----------------------------------------

%parametros - constantes
df      = fs/N;                         % incremento de frecuencia. Escala lineal
Nmax    = N/2;                          % max fft indice
fmax    = fs/2;                         % max frec.
melmax  = freq2mel(fmax);               % max mel frec.
melinc  = melmax / (nofChannels + 1);   % increment frec. en escala de mel

%frecuencias centrales. Mel
melcenters = (1:nofChannels) .* melinc;

%frecuencias centrales. Hz
fcenters = mel2freq(melcenters);

%calcula los BW de los filtros. Escala lineal
startfreq   = [0 , fcenters(1:(nofChannels-1))];
endfreq     = [fcenters(2:nofChannels) , fmax];
bandwidth   = endfreq - startfreq ;

%quantize into FFT indices
indexcenter = round(fcenters ./df);

%calcular frecuencias resultantes.compute resulting frequencies
fftfreq = indexcenter.*df;

%calcular el punto inicial de los indices de las ventanas.
indexstart = [1 , indexcenter(1:nofChannels-1)];

%compute stop indices of windows
indexstop = [indexcenter(2:nofChannels),Nmax];

%calcular bw (numero de indice por ventana)
idxbw = (indexstop - indexstart)+1;

%calcular BW [Hz]
FFTbandwidth = idxbw.*df;

%calcular el quantization error resulatante
diff = fcenters - fftfreq;


%calcula la matriz de coeficientes de los filtro triangulares
W = zeros(nofChannels,Nmax);
for c = 1:nofChannels
    %left ramp
    increment = 1.0/(indexcenter(c) - indexstart(c));
    for i = indexstart(c):indexcenter(c)
        W(c,i) = (i - indexstart(c))*increment;
    end
    %right ramp
    decrement = 1.0/(indexstop(c) - indexcenter(c));
    for i = indexcenter(c):indexstop(c)
       W(c,i) = 1.0 - ((i - indexcenter(c))*decrement);
    end
end

%Ws=sparse(W); %guarda la matriz sin ceros
%save matrix
%save 'melmatrix' W;

%plot results

%plot filter frequencies

% subplot(2,3,1);
% figure;
% plot(fcenters,melcenters,'-o');
% title('Frecuencias de Filtros en escala de Mel');
% xlabel('f[Hz]');
% ylabel('f[mel]');
% grid on;
%  
% %plot channel numbers and center frequency FFT indices
% %figure;
% subplot(2,3,2);
% plot(0:nofChannels-1,indexcenter,'-o');
% title('FFT indices');
% axis([0,nofChannels-1,-Inf,Inf]);
% xlabel('Canal del filtro');
% ylabel('FFT indice');
% grid on;
% 
% %plot frequency quantization error
% figure;
% subplot(2,3,3);
% plot(0:nofChannels-1,(diff),'-o');
% title('Frequency quantization error');
% axis([0,nofChannels-1,-Inf,Inf]);
% xlabel('Canal del filtro');
% ylabel('f_{centro} - f_{fft} [Hz]');
% grid on;
% 
% 
% %plot filter bandwidth
% subplot(2,3,4)
% plot(0:nofChannels-1,bandwidth,'-o');
% axis([0,nofChannels-1,-Inf,Inf]);
% title('Ancho de Banda del filtro');
% xlabel('Canal del filtro');
% ylabel('Ancho de Banda [Hz]');
% grid on;
% 
% %plot indexbandwidth
% subplot(2,3,5);
% plot(0:nofChannels-1,idxbw,'-o');
% title('Ancho de Banda del filtro FFT[indices]');
% axis([0,nofChannels-1,-Inf,Inf]);
% xlabel('Canal del filtro');
% ylabel('N° de FFT indices');
% grid on;
% 
% %plot bandwidth error
% subplot(2,3,6);
% plot(0:nofChannels-1,(bandwidth - FFTbandwidth),'-o');
% title('Bandwidth quantization error');
% axis([0,nofChannels-1,-Inf,Inf]);
% xlabel('Canal del filtro', 'FontSize',12);
% ylabel('Error de Ancho de Banda (BW error)[Hz]', 'FontSize',12);
% grid on;
% 
% %plot filter coefficients W
% figure;
% imagesc(1-W);
% grid on;
% title(['Matriz de coeficientes de los filtros Mel. N:',sprintf('%s',num2str(nofChannels))], 'FontSize',12);
% xlabel('Indice de FFT', 'FontSize',12);
% ylabel('Canal del filtro','FontSize',12');
% colormap(gray);
% 
% figure;
% subplot(2,1,1)
% hold on
% for f=1:nofChannels
%         plot(W(f,:),'LineWidth',2, 'Color', [0.04,0.52,0.78]);
% end
% 
% title(['Banco de ' sprintf('%s',num2str(nofChannels)) ' filtros Mel'],'FontSize',12);
% 
% subplot(2,1,2)
% stems = zeros(size(fcenters));
% plot(fcenters,stems,'Marker', 'o', 'MarkerSize',8, 'MarkerFaceColor', [0 0 1])
% title('Frecuencia Central [Hz]', 'FontSize',12)